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►^ Abstract 

>0 In this paper we compute two important inforraation-theorctic quantities which arise in 

~V the application of multiple-input multiple-output (MIMO) antenna wireless communication 

^--^ systems: the distribution of the mutual information of multi-antenna Gaussian channels, and 

^* the Gallager random coding upper bound on the error probability achievable by finite-length 

^-^ channel codes. It turns out that the mathematical problem underpinning both quantities is 

f~^ the computation of certain Hankel determinants generated by deformed versions of classical 

T-H weight functions. For single-user MIMO systems, it is a deformed Laguerre weight, whereas 

L| for multi-user MIMO systems it is a deformed Jacobi weight. We apply two different methods 

• ^H to characterize each of these Hankel determinants. First, we employ the ladder operators 

/\^ of the corresponding monic orthogonal polynomials to give an exact characterization of the 

Hankel determinants in terms of Painleve differential equations. This turns out to be a 
Painlcvc V for the single-user MIMO scenario and a Painleve VI for the multi user scenario. 
We then employ Coulomb fluid methods to derive new closed-form approximations for the 
Hankel determinants which, although formally valid for large matrix dimensions, are shown 
to give accurate results for both the MIMO mutual information distribution and the error 
exponent even when the matrix dimensions are small. Focusing on the single-user mutual 
information distribution, we then employ both the exact Painleve representation and the 
Coulomb fluid approximation to yield deeper insights into the scaling behavior in terms of 
the number of antennas and signal-to-noise ratio. Among other things, these results allow us 
to study the asymptotic Gaussianity of the distribution as the number of antennas increase, 
and to explicitly compute the correction terms to the mean, variance, and higher order 
cumulants. 



1 Introduction and Preliminaries 

Over the past decade, multiple-input multiple-output (MIMO) systems have been at the forefront 
of wireless communications research and development. Such systems, which employ multiple 
antennas at both the transmitter and receiver, have revolutionized the wireless industry and now 
form the basis of most emerging wireless standards, such as next-generation wireless local area 
networks (WLAN) and cellular mobile networks. The main reason for this explosion of interest 
is due to the independent discoveries of Telatar [1] and Foschini [2], which demonstrated that the 
fundamental information-theoretic capacity of MIMO systems grows linearly with the number of 
antennas. This is important, since traditional methods for increasing capacity, which typically 
increase valuable system resources such as the transmission power, yield only a logarithmic 
capacity increase. Indeed, MIMO is now widely recognized as a key technology for meeting the 
ever-increasing demands for higher-rate data-oriented wireless communications applications and 
services. 

There are various metrics for characterizing the fundamental information-theoretic limits of 
MIMO systems, each of which have relevance depending on the specific wireless communica- 
tion configuration. The most commonly studied metric is the so-called ergodic capacity, which 
specifies the maximum achievable average mutual information between the transmitter and re- 
ceiver. This measure, which assumes that there is sufficient dynamics in the system such that a 
user's codeword may span a large number of "independent channels", has been studied exten- 
sively over the past decade (see e.g., [1-10] and the tutorial discussion [11]). A closely related 
measure which has received far less attention is the outage capacity. This measure is impor- 
tant for characterizing the communication limits of systems which are not highly dynamic (for 
example, WLANs), or systems which have stringent delay requirements. Compared with the 
ergodic capacity, the analysis of outage capacity is much more complicated since it requires the 
distribution of the mutual information between the transmitter and receiver, rather than simply 
the average value. Another important information-theoretic quantity of practical interest is the 
error probability achievable with block-coding schemes of a given length and rate. Whilst an 



exact characterization of this quantity is not tractable in general, methods have been proposed 
by Gallager for upper bounding this quantity [12]. Such bounds have been well-established for 
single- antenna systems, however much less is known for MIMO. 

In this paper, we aim to present new methods for studying the outage capacity and the error 
probability of MIMO systems. As we will see, in both cases the fundamental mathematical 
problem of interest boils down to characterizing the Hankel determinant 



D„ = det (^i+j)" (1.1) 



generated from the moments of a certain weight function w{x) 



lik-= x^w{x)dx, A; = 0,1,2,... (1.2) 

with J denoting the support of the weight. The specific form of the weight depends on the 
MIMO configuration under investigation. We will consider two important MIMO configurations 
pertaining to single-user and multi-user MIMO systems respectively. In the single-user case, the 
weight function is shown to be a deformation of the classical Laguerre weight given by 

u;(x) = x" e"'^ (x -h t)^, 0<x<oo, a>0, t>0, (1.3) 

whereas in the multi-user case the weight function is found to be a deformation of the (shifted) 
classical Jacobi weight given by 

w;(x) = x"i(l-a;)"' (^^) , xG[0,l], ai > 0, as > 0, t > 0. (1.4) 

We employ two different methods from random matrix theory for evaluating the correspond- 
ing Hankel determinants, both of which have been used to great effect by the mathematical 
physics community. First, we derive exact expressions for the Hankel determinants by employ- 
ing the theory of orthogonal polynomials and their corresponding ladder operators. There exists 
extensive literature on this subject; for example [13-20]. See also [21-26] for recent applications 



of ladder operators to Hermitian matrix ensembles. 

Despite these prior contributions, this approach has not been employed to problems in 
wireless communications and information theory. Through the ladder operator framework, we 
demonstrate that the Hankel determinant generated from the deformed Laguerre weight (1.3) 
has a simple representation involving the Painleve V differential equation. We also derive a dis- 
crete difference equation representation, which allows the Hankel determinant to be numerically 
evaluated iteratively, in terms of the Hankel matrix dimension. For the deformed Jacobi weight( 
1.4), upon establishing the connection to the multi-user MIMO system model of interest, we 
employ a known result derived using the ladder operator method in [27] to represent the corre- 
sponding Hankel determinant in terms of the Painleve VI differential equation. This presents 
a new connection between the Painleve VI differential equation and the information theory of 
multi-user MIMO communication systems. 

In addition to deriving exact expressions for the Hankel determinants, we also derive ap- 
proximations for these determinants by employing the general linear statistics results from [28], 
derived based on Dyson's Coulomb fluid interpretation [29-31]. These asymptotic results are 
essentially the Hankel analog of Szego's strong limit theorem on the asymptotic characterization 
of large Toeplitz determinants, a component of which appeared as early as 1919 [32], where Han- 
kel determinants generated by compactly supported weights were studied. We refer the reader 
to [33-35] for more related contributions. 

As with the ladder operator approach. Coulomb fluid techniques have been used very suc- 
cessfully in the context of mathematical physics, however, they have yet to be adopted by the 
wireless communications and information theory community. (An exception is the very recent 
submitted work [36], which employed the Coulomb fluid framework but not the linear statistics 
method to characterize the large-antenna behavior of the single- user MIMO channel capacity.) 
A key advantage of the Coulomb fluid methodology is that the expressions which are obtained 
are in closed-form, and are simpler than those obtained via exact methods. Moreover, although 
formally valid for large matrix dimensions, we find that the approximations are remarkably 
accurate for even very small matrix dimensions (e.g., 2x2). By adopting this Coulomb fluid 



framework, our results also establish the Gaussian behavior of the channel capacity as the num- 
ber of antennas grow large, for both the single-user and multi-user MIMO systems. This point 
has been derived previously for the single- user MIMO case, using different methods from [37], 
however we do not believe that it has been established for the multi-user MIMO scenario. 

After deriving the exact Painleve and approximate Coulomb fluid representations, we then 
employ both results to gain useful insights into the mutual information distribution. In particu- 
lar, focusing on the single-user MIMO scenario (i.e., deformed Laguerre weight) with n transmit 
and receive antennas, we demonstrate that the Coulomb fluid approximation in fact gives the 
exact distribution of the mutual information, corresponding to a Gaussian, to leading order in 
n. We also use the Painleve V equation to compute the large-n correction terms for the mean, 
variance, and third cumulant. Among other things, these results demonstrate the interesting 
phenomenon that the distribution deviates from Gaussian as the signal-to-noise ratio (SNR) 
increases. The sensitivities of the mean, variance, and third moment, with respect to the SNR 
are also examined. 

The paper is organized as follows. In the rest of this section, we present a detailed discussion 
of the linear models which underpin the single-user and multi-user MIMO wireless communica- 
tion scenarios of interest. We also introduce the fundamental information-theoretic measures of 
outage capacity and error probability, and establish important connections with Hankel determi- 
nants generated from deformed Laguerre and Jacobi weights. Then, in Section 2, we introduce 
the ladder operators of orthogonal polynomials and their associated compatibility conditions, 
which provide the key ingredients for establishing an exact finite n characterization of the Hankel 
determinants in terms of Painleve differential equations. In Section 3, we introduce, for large 
n, the probability density of a class of random variables called linear statistics. These results, 
based on Dyson's Coulomb fluid interpretation [29] and developed further in [28,30,31,38], are 
very general and embrace a wide class of random matrix models. By employing these general 
results, we derive closed-form approximations for the mutual information distribution for the 
Hankel determinants generated by both the deformed Laguerre and Jacobi weights. These re- 
sults permit very fast evaluation of the error probabilities of single-user and multi-user MIMO 



systems, whilst also demonstrating the Gaussian behavior of the mutual information distribution 
for large n. In Section 4, focusing on the single-user MIMO scenario, the Coulomb results are 
compared, for large n, with the solutions of the continuous cr-form of the Painleve V. We also 
compute closed-form asymptotic expressions for the recurrence coefficients corresponding to the 
deformed Laguerre weight, which are basic variables in our theory. Our analysis will involve a 
number of complicated integrals, which are derived and tabulated in the Appendix. 

1.1 Information Theory of MIMO Wireless Systems 

Consider a MIMO communication system with rit transmit and n^ receive antennas. The linear 
model relating the input (transmitted) signal vector x„jxi G C"' and output (received) signal 
vector Yrirxi £ C"'' takes the form 

y = Hx + n. (1.5) 

Here, n„^xi £ C"'' is a complex Gaussian vector with zero mean and covariance £'(nn"l") = Q„. 
This covariance matrix can account for the effects of both receiver noise as well as multi-user 
interference, and as such, the selection of Q„ will distinguish between the single-user and multi- 
user MIMO models which we consider subsequently. The matrix H G £nrxnt ^ referred to as 
the channel matrix^ represents the wireless fading coefficients between each transmit and receive 
antenna. This matrix is assumed to be known to the receiver^, but not to the transmitter. 
The channel is modeled stochastically, with distribution depending on the specific wireless en- 
vironment. Under the realistic assumption that there are sufficient scatterers surrounding the 
transmit and receive terminals, the channel matrix H is accurately modeled according to a com- 
plex Gaussian distribution with independent and identically distributed (i.i.d.) elements having 
zero mean and unit variance. The transmitted signal x is designed to meet a power constraint: 



^(xfx) < P . (1.6) 



^In practice, this information can be obtained using standard channel estimation techniques. 



1.1.1 Ergodic and Outage Capacity 

The Shannon capacity governs the ultimate hmits of communication systems. More specifically, 
this measure defines the highest data rate that can be achieved with negligible errors by any 
transmission scheme. As such, the Shannon capacity forms a benchmark for the design of 
practical transmission technologies. Mathematically, the Shannon capacity is defined in terms 
of the mutual information between the input and output signals, which for the MIMO linear 
model (1.5) is given by: 

/(x;y|H) = n{y\ll) - n{y\^,U) 

= niy\ll)-n{n) (1.7) 

with 'H(y|H) denoting the conditional entropy of y, defined in terms of its density p(y|H) as 
follows: 

^(y|H) = E{- logp) :=- f p(y|H) logp(y|H)dy. (1.8) 

In general, there are two important measures of capacity — the ergodic capacity and the outage 
capacity. The ergodic capacity is the relevant measure for applications with highly dynamic 
channels; for example, high-mobility wireless applications, in which case the channel matrix 
H varies quickly over time, and therefore each transmission codeword sees a large number of 
"independent" channel realizations. The ergodic capacity is defined as: 

C = max£;H(/(x;y|H)) (1.9) 

p(x) 

where the maximum is taken over all densities p(x) of the input vector x, subject to the constraint 
(1.6). It has been proven in [1] that the optimal input density p*(x) is multi-variate complex 
Gaussian with zero mean. As such, the mutual information /(x;y|H) becomes 

I(x;y|H) = logdet (v + HQHtQ^i) (1.10) 



where Q = £'(xx') is the input signal covariance. The capacity (1.9) can therefore be reposed 
as an optimization over Q as 

C = max£;H(/(x;y|H)) (1.11) 

subject to tr(Q) < P. This quantity has been studied extensively over the past decade. Indeed, 
methods have been proposed to calculate the optimal Q for a range of MIMO channel models of 
interest [39, 40], and the expectation has been characterized through the use of both large- and 
finite-dimensional random matrix theory [1-6,9, 11,41]. For the case which we consider, where 
H is assumed to be complex Gaussian with independent zero- mean unit- variance entries, the 
capacity-achieving input covariance has been derived in [1] and is given by 

Q* = -In, . (1.12) 

nt 

In practical terms, this implies that the capacity is achieved by sending independent Gaussian 
signals from each of the transmit antennas with equal power. 

Compared with the ergodic capacity, the outage capacity is a much more difficult problem 
and has received far less attention. In contrast to the ergodic capacity, this measure is suitable 
for wireless applications with low mobility (e.g., wireless local area networks), in which case 
the channel is assumed fixed during the transmission of a codeword. The outage capacity 
corresponding to an outage probability Pout is defined as the transmission rate which can be 
supported by (1 — Pout) x 100% of the channel realizations^. Although the optimal distribution of 
x which maximizes this quantity is unknown, a sensible choice is to employ the same distribution 
as that which achieves the ergodic capacity, i.e., p*(x). In this case, the outage capacity Cout 
satisfies: 

Pout(Cout) = Pr(/(x;y)<Cout) 



By supported, we mean that the mutual information for a given channel realization is greater than the 
transmission rate. 



= Pr Aogdetfl„, + ^HHtQ;i') <Cout') (1.13) 

with Q* denoting the input covariance which maximizes (1.11). Compared with (1.11), the 
quantity on the right-hand side (r.h.s.) of (1.13) is much more difficult to characterize since 
it involves the entire distribution of the logdet(-) random variable, rather than simply the 
expected value. As such, a common approach has been to compute the first few moments of the 
distribution, and then use these to obtain a Gaussian approximation (see e.g., [4,6]). 
The outage probability (1.13) can be calculated via 

Pout(Cout) = — / Mijoj) ; du (1.14) 

where Ai{-) denotes the moment generating function of the mutual information, taking the form 

M{X) := ^H(exp(A/(x;y|H))) 

= ^HldetTv + ^HHtQ^A J, (1.15) 

and j := \/— T. 

Alternatively, the Chernoff bound may be employed, as in [9], to give 

Pout(aut) < M{-s)e'^°-^ (1.16) 

for any s > 0. Let 

Gis):=logMi-s), 

and rewrite the moment generating function in an obvious alternative form: 



Then a straightforward computation shows that 

^ f e-^'i(^) [q{X)]^f^{dX) ( f e--^Wg(X)^(dX) \ ' 
^'' Je-'i(^)lJ.{dX) \ Je-'i(^)n{dX) J 

If we interpret 

e-^'?(^V(rf^) 

as a probability measure, then for all real s, 



where 



F{X): 



G"{s) = {q-qy>0, 



_ jF{X)e-'^(^')fi{dX) 



Je-^i(^)fi{dX) ' 
Since G{s) is convex we can minimize the Chernoff bound by minimizing 

Gis) + s Gout 

and minimum of the above is 

1 {s* G (s*) - Cout) 1 *2W// *x 

"2 GVO 2" ^^'^' 

where G{s) is minimized at s* and we have assumed that s* > 0. 

1.1.2 Upper Bound on Error Probability 

Whilst the capacity is an important fundamental quantity, it is difficult to approach in practice 
since it requires infinitely long codewords and receivers with unbounded complexity. Thus, an 
important question is to determine what data rates R are achievable for practical channel coding 
strategies with fixed and finite codeword length N, subject to a requirement on the tolerated 

10 



error probability Pc- One method for addressing this problem was proposed by Gallager [12], 
where a general upper bound on P^ was derived, assuming that randomly- selected length-A^ 
block codes of rate R were employed, along with maximum-likelihood receivers. For the MIMO 
transmission model (1.5), assuming that the channel is memoryless (i.e., each transmission 
period sees an independent realization of H), the Gallager random coding upper bound on error 
probability is expressed as [1] 



Pc(iV,i?)<e-^^'-(«), 



(1.17) 



where Er{R) is referred to as the error exponent, which is independent of N. This function is 
given by 



EriR) = maxo<p<i {Eo{p) - pR}, 



(1.18) 



with 



Eo{p) = -log Eu 



C"t 



p(x){p(y|x,H)}i/(i+^)dx 



(i+p) 



dy 



where 



p(y|x,H) 



exp[ 



Hx)tQ-i(y 



Hx)] 



det(7rQr 



(1.19) 



(1.20) 



and p(x) denotes the density of the input signal x. To maximize the error exponent and thus 
minimize the error probability, p(x) should be selected so as to maximize Eq{p). Evaluating 
this optimal input distribution is very challenging, and a sensible (and more tractable) choice 
is to assume that x takes the capacity-achieving distribution p*(x) presented in the previous 
section; i.e., it is zero-mean complex Gaussian with covariance Q*. In fact, it was shown in [1] 
that under the assumption that x is Gaussian, p*(x) is optimum in terms of maximizing Eq{p). 
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With this input distribution, it is easy to see that (1.19) particularizes to: 



Eo{p) 



log E^ 



det I„, + 



HHtQ^i 



nt{l+p) 



(1.21) 



We remark that whilst more refined bounds on the error probability compared with (1.17) have 
also been derived [12] , these more elaborate bounds still yield the same underlying mathematical 
problem as that posed in (1.21). 

1.2 Single-User MIMO and the Deformed Laguerre Weight 

Single-user MIMO systems embrace a class of coordinated wireless networks for which all trans- 
missions are scheduled in an orthogonal manner (e.g., in orthogonal time-slots or orthogonal 
frequency bands), and as such, transmissions do not interfere with one another. In this scenario, 
the transmitted and received signals conform to the linear model (1.5), with n simply refiecting 
the receiver noise. This noise is spatially uncorrelated (across antennas), and without loss of 
generality it has covariance 



Q. 



(1.22) 



Due to the normalization of the trace of Q„, the transmit power P also represents the SNR. 

The key quantity of interest for the outage capacity is the moment generating function, which 
in this case particularizes to 



M{\) = Eu 



det 



In^ + -HHt 



(1.23) 



with 



t:-- 



P 



(1.24) 



Similarly, the error exponent (1.19) admits the same form, but with substitutions A = —p and 
t = {l + p)nt/ P. In the following, we will focus our discussion on the moment generating function 
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(1.23), keeping in mind that the apphcation to the error probabihty is immediate. 
Let 

m := max{nr,nt}, n := m\n{nr , nt} , a := m — n 

and define 



W 



HH^', llr < fit 

H''^H, Ur > nt 



The matrix W is a complex Wishart random matrix with positive eigenvalues denoted by 
{xi}f^i. It is well known that the joint probability density function of the eigenvalues read 



p(xi,X2,...,x„) oc JJ-u;Lag(a:i) JJ (xj-Xkf, 



i=l 



i<j<k<n 



(1.25) 



where Xj G [0, oo), and tyLag(") is the classical Laguerre weight 



w'Lag(x) = x"e ^. 



With these definitions, using the identity 



det(I + AB) = det(I + BA), 



(1.26) 



we can evaluate the moment generating function as 



M{X) = E 



det ( I„ + -W 



n(i+ 

LA:=1 



¥)■ 



'DLl (1+ ^) ^i'Lag(xfc)dXfc 



Iwi Y\i<Mi - Xjf Uk=l W'L^g{xk)dXk 



(1.27) 
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The remaining integrals are evaluated in determinant form via the Andreief-Heine identity: 

If " 

det{ni+j)1j^Q = Dn[w] = — TT {xj - Xif TT w{xk)dxk, (1.28) 



{a,b)-~ l<i<j<n fc=l 

where 



Hi:= I x'w{x)dx, i = 0,1,2,... (1.29) 



are moments of the weight w. See [42] for a recent exposition on this and other related matters. 
Obviously the moments would depend on the parameters which may appear in the weight. 
With this identity, we immediately obtain 

where 

Dn{t, A) = det {fi,+jit, A))^T^o (1.31) 

is the Hankel determinant generated from the deformed Laguerre weight 

w'dLag(a;) = w'dLag(a;, t. A) := (x + t)^WLi,gix), t>0 (1.32) 

with moments 

/>oo 

;Ufc(t, A):=/ x''wdi^a,g{x)dx , A; = 0, 1, 2, . . . . (1.33) 

Jo 

Remark 1 The factor D„(t,0) is simply the Hankel determinant generated from the non- 
deformed Laguerre weight, wi,3,g{x), which can be computed exactly in terms of the Barnes G— 
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function as 

i^n(M)=i^„,.Ka,]= ^^"+^^^: + " + ^\ G(l) = l. (1.34) 

Similarly, for t = and fixed A, it follows that 

Dn{0, A) = Dn,a+x[wLa.g\ ■ 

Remark 2 The moments are expressed in terms of the Kummer function of the second kind 
U{a; b; z) as follows: 

/ife(t,A) := / x^+'^Cx + 1)^ e-^dx 

= t"+^+'=+i Y{a + k + 1) U{a + fc + 1; q + A + fe + 2; t), /c = 0, 1, 2, ... (1.35) 

Note that this expression has been previously reported in [43], and an alternative representation 
given in [44]. Whilst this identity, combined with (1.34) and (1.30), gives a "closed-form" 
determinantal representation for the moment generating function, it does not provide useful 
insights and it also becomes unwieldy to evaluate if the number of antennas become large. 
Moreover, in this form, it is not amenable to further analysis, such as determining the asymptotic 
scaling properties. To overcome these problems, in Section 2 we employ the theory of orthogonal 
polynomials and their corresponding ladder operators to provide a more useful characterization, 
where we express the Hankel determinant generated from the deformed Laguerre weight in terms 
of the classical Painleve V differential equation. 

An alternative characterization for the moment generating function which will also prove 
useful is derived as follows. Starting with (1.27), and applying the transformations Xj — )■ nxj, 
i = 1, . . . , n, we obtain 

JM"n,<,(^.-^.f nLi {T+ x,f x^e— Mxfc 
M{X) = T~"^ + .' -:^,= (1.36) 
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where T := /3i/-P, with /3i = nt/n. Equivalently, 



-nA^«(A) 



where 



•^(^' = ^ iM "■'" 



c ( '^ \ " 

Zn{X)= exp -$(xi,...,x„) + Aj];iog(T + x,) m^^fc (^-3^) 

•^'*+ \ j=i / fc=i 



with 



$(xi,...,x„) = -2 ^ log|xj -Xfcl +n^(x£- /31ogX£) (1.39) 

^<j<k<n 1=1 

and we have defined /3 := m/n — 1. 

This representation will be important for deriving an approximation for the moment gener- 
ating function based on the Coulomb fluid approach in Section 3. 

1.3 Multi-User MIMO and the Deformed Jacobi Weight 

In this section we introduce the second communication scenario which we will consider, corre- 
sponding to a multi-user MIMO system. Such systems are highly relevant for cellular mobile 
applications, where the transmit-receive communication channel is impaired by interference from 
other users operating within the same frequency bandwidth. In fact, the key issue of interference 
presents one of the most important challenges in the successful deployment of MIMO in practice 
[45]. We will focus on the practical interference-limited scenario, where the receiver noise is neg- 
ligible compared with the multi-user interference. We assume that there are K interferers, each 
equipped with ut antennas, and transmitting (capacity achieving) independent Gaussian signals 
with power Pj/rit out of each transmit antenna. We make the common assumption (see e.g., 
[46,47]), that the interferers have equal-power, which is valid when the interferers are located 
at similar distances to the receiver. Moreover, as discussed in [46], if this assumption is not 
met, then the equal power assumption leads to a lower bound on performance. We also assume 
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that rir < Tit, which is appropriate for modehng the downhnk communication (base-station to 
mobile transmission) of a interference-hmited cellular system, where the number of transmit an- 
tennas deployed by the base-station may be large, but the number of receive antennas is highly 
restricted due to limited space constraints of the cellular mobile device. 

Under the assumptions discussed above, the noise term n in (1.5) takes the form 

K 



n = J^H,Xi (1.40) 



where Hj and Xj denote the channel matrix and input vector for the ith interferer respectively. 
These are assumed independent across i, and independent of H and x. It is convenient to write 
this in the stacked form 

n = H/x/ (1.41) 



where 



Hi = [Hi, • • • , Hi^] , X7 = [xi, • • • , jckV . (1.42) 

Here, H/ G C"''^^"' is complex Gaussian with independent zero- mean unit- variance entries, 
whilst x/ € C * is complex Gaussian with independent zero-mean entries having variance 
P/nt- The noise covariance matrix, conditioned on H/, is therefore given by 

Ctn = — H/HI . (1.43) 

lit 

In this case, the moment generating function of the capacity particularizes to 

7W(A) = ^H det Tv + -^HHt(H7H|)-i') . (1.44) 

Similarly, the error exponent (1-19) admits the same form, but with the substitution A = —p and 
replacing P/ with (1 + p)Pi. Here, HH^ and H/H| are independent complex central Wishart. 
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As for the deformed Laguerre case, in the fohowing we wih focus our discussion on the moment 
generating function (1.44), keeping in mind that the apphcation to the error probabihty is 
immediate. 

For consistency with previous notation, let us define: 

n := rir, mi := rif, 1712 '■= Krit, ai := mi — n, 02 '■= 1^2 — n ■ (1-45) 

Generahzing [48, pp. 312-314] from real to complex matrices, we find that the joint probability 
density function of the eigenvalues /i, ..., /„ of the random matrix HHT(H/Hj)~"'^ is given by 

P{fl, /2, ...,/„) OC JJ Arm+m2 H Ui-fjf, 

k=l ^ ■' ' l<i<j<n 

where fk G (0, 00), k = 1,2, ..., n. With the change of variables^ 

f Xk 

Jk 



1-Xk 

the above density becomes 

n 
/(X1,X2,...,X„) DC JJu;jac(Xfc) JJ {Xi-Xjf, (1.46) 

fc=i i<*<i<™ 

where x^ G (0, 1), k = 1,2, ...,n, and tyjac(0 denotes the shifted classical Jacobi weight 

WiUx) = x'''{l-xr\ xG[0,l]. (1.47) 

Note that the classical Jacobi weight has the form 

*Jac(x) = (l-x)"i(l+^)"', xG[-l,l] (1-48) 



'We would like to thank Iain Johnstone for pointing this out. 
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and therefore 



^Jac(x) = ^^-l^V J^\ XG[0,1]. 



(1.49) 



The moment generating function can be evaluated as 



Mix) = E 



E 



.k=l 



P 



n 1+^/^ 



Pi' 



.fc=l ^ 



Pjl- Xk 



t~'^^E 



n 

.fc=i ^ 



t + Xk 



Xk 



L—nX 



/(0,1)" Ul<i<j<niXj - Xi? riLl ^Jac(Xfc) (^^ j dx, 
J(0,1)" lll<i<j<n(^i ~ ^i) [[k=l^i!^c\^k)dXk 



where 



t 



Pi 



P-Pi 

We note here that the t variable has no n dependence when expressed in terms of Pj/{P—Pj), 
rather unlike the single user MIMO case. This is mathematically reasonable as the Jacobi weight 
is compactly supported. 

Applying the Andreief-Heine identity (1.28), we obtain 



M{X) = t 



-nxDn{t,\,mi,m2) 
L>„(t,0,mi,m2) 



(1.50) 



where 



Dn{t,X,mi,m2) = det {fii+j{t, X, 1712, m2))'^j^Q 



(1.51) 



is the Hankel determinant generated from the deformed (shifted) Jacobi weight 



Wdj!,c{x) = Wdji^cix, t, A) := wjs.cix) 



x + t 
1 — X 



(1.52) 
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with moments 

^fc(i, A,mi,m2) := / x''wdj!,c{x)dx , A; = 0,1,2,... (1.53) 

Jo 

Remark 3 The factor Dn{t,0,mi,m2) is simply the Hankel determinant generated from the 
(shifted) non-deformed Jacobi weight, wjacix), which can be computed exactly in terms of the 
Barnes G— function. To this end, we apply the transformation (1.49) to give 



1 1 



„ n 

■ ■J(-'-^i-) l<i<j<n k=l 

(27r)'^ r (£l+|2+l) (-2 (^01+02+1^ ^2 (^01+02 ^ ^^ 

4"("i+"2+n) Q (q,^ + q!2 + 1) G (ai + 1) G (a2 + 1) 
G{n + l)G{n + ai + l)G{n + 02 + l)G{n + ai + 02 + 1) 
Q2 (^ ^ gi+ag+l^ Q2 (^ ^ «l±«2 ^ 1-) r (n + "i+"+2 +iy 



(1.54) 



where, to obtain the last equality, we have invoked [35, Eq. (1.6)]. 
Similarly, for t = and fixed A, it follows that 

^^(O, A, mi, 7712) = Dn,ai+\,a2-x[wji,c] ■ (1-55) 

Remark 4 The moments can be expressed in terms of the hypergeometric function 2-^1 (•) as 
follows: 

Hk{t, A, 1711,1712) = t^T{l + k + ai)r(l + Q2 - A) 2i^i(l + fc + Qi, -A, 2 + fc + ai + a2 - A, -1/t). 

This relation was also pointed out in [49]. Whilst this identity, combined with (1.55) and (1.50), 
gives a "closed-form" determinantal representation for the moment generating function, it does 
not provide useful insights and it also becomes unwieldy to evaluate if the number of antennas 
become large. To overcome these problems, in Section 2 we employ the theory of orthogonal 
polynomials and their corresponding ladder operators to provide a more useful characterization. 
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where we express the Hankel determinant generated from the deformed Jacobi weight in terms 
of the classical Painleve VI differential equation. 

Similar to (1.37), it will also be useful to note the following equivalent representation for the 
moment generating function, 



where T = t, 



[ exp(-$(xi,...,x„) + AVlogf^^^H "TTdxfc (1.57) 

J (0,1)- \ ti \^-^Jjt=\ 



and 



n 



$(xi,...,x„) = -2 ^ log\xi - Xj\-n'^ {ipilogXi + ip2log{l- xe)) (1.58) 



with 



ipi = 1, ^2 = 1 • (1-59) 

n n 

Once again, this representation will be critical for employing the Coulomb fluid methodology in 
Section 3. 



2 Exact Characterization via the Ladder Operator Method 

In the previous section, we demonstrated strong relationships between two information measures 
of MIMO channels, namely, the outage capacity and the error probability, and certain Hankel 
determinants. For single-user MIMO systems, the Hankel determinant of interest was generated 
via the moments of a deformation of the Laguerre weight, whereas for multi-user MIMO systems, 
it was generated via the moments of a deformed Jacobi weight. In this section, we present an 
exact characterization of these Hankel determinants by employing the theory of orthogonal 
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polynomials and their corresponding "raising and lowering" ladder operators. Before presenting 
the main results, we first introduce some preliminary material which will prove useful. 

2.1 Preliminaries of Orthogonal Polynomials and their Ladder Operators 

Here we provide a brief discussion to highlight the connections between orthogonal polynomials 
and the Hankel determinants of interest, as well as provide basic properties of orthogonal poly- 
nomials which will be needed. From (1.25) and (1.46), the joint eigenvalue distributions arising 
in the single-user and multi-user MIMO scenarios admit the generic form 

n 

P{yi, • • • , yn) oc Jl w{yk) Y{ ^yj ~ ^»)^' 2/j G {a, b) 

fc=l l<j<j<n 

with w{-) denoting the weight. Moreover, in both cases, the key quantity of interest takes the 
generic form 

If "^ f f^ \" 

^' Aa,ft)" l<y<n Ll VA y„=0 

where 

w{y) = w{y)g{y) 

denotes the deformed weight. As a direct consequence of the Vandermonde determinant factor, 
by applying elementary operations, this Hankel determinant can be equivalently expressed as 

Dn = Aei( I P^{y)P,{y)w{y)dy\ 

where Pj{-) represents any monic polynomial of degree j, written as 

P,{z) = z^ + p,{j)z^~' + ... (2.1) 
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If we orthogonalize the polynomial sequence {Pn{y)} with respect to w{y) over the interval [a, 



I.e., 



b 

Pi{y)Pj{y)w{y)dy = hA,o: hj = 0, 1,2, ... (2.2) 



with hi denoting the square of the L? norm of Pi , then the Hankel determinant D„ evaluates to 

Dn=Whk . 
k=0 

Thus, we clearly see that the problem of computing the Hankel determinants of interest becomes 
one of characterizing the class of polynomials which are orthogonal with respect to the deformed 
weight function. To attack this problem, we require some definitions and tools, as given below. 
We start by noting that if 

|Ui := / y^w{y)dy, 
J a 

exists for all i = 0,1,2,..., then the theory of orthogonal polynomials states that Pn{y) for 
n = 0, 1, 2, ... satisfies the three term recurrence relations, 

zPn{z) = Pn+l{z) + an Pn{z) + /3„ Pn-l{z). (2.3) 

The above sequence of polynomials can be generated from the orthogonality conditions, the 
recurrence relations and the initial conditions, 

Po{z) = l, /3oP-i(z) = 0. 

For example 

Fi[z) = z — aQ = z . 

^0 
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Substituting (2.1) into the recurrence relations, an easy computation shows that 

Pi(n)-pi(n + l) = a„, (2.4) 

with Pi(0) := 0. A telescopic sum of (2.4) gives 

an = -Y.PiiJ). (2.5) 

i=o 

From the recurrence relation (2.3) and the orthogonality relations (2.2), we find 

/3n=T^. (2.6) 

We shall see that Pi{n) plays an important role in later developments. For more information 
on orthogonal polynomials, we give reference to Szego's treatise [50]. 

Next, we present three lemmas, which represent the ladder operators of orthogonal polyno- 
mials, as well as some supplementary conditions. Note that these results have been known for 
quite sometime; we reproduce them here for the convenience of the reader using the notation 
of [24], where one can also find a list of references to the literature. We would like to mention 
here that Magnus [51] was perhaps the first to apply these lemmas — albeit in a slightly different 
form-to random matrix theory and the derivation of Painleve equations. It should be mentioned 
that Tracy and Widom also made use of the compatibility conditions in their systematic study 
of finite n matrix models [52]. See also [53, 54]. In presenting these lemmas, we use the following 
definition: 

V = —logw. 

Lemma 1 Suppose v = — log w has a derivative in some Lipshitz class with positive exponent. 
The lowering and raising operators satisfy the following: 

P'^{Z) = -Bn{z)Pn{z)+PnAn{z)Pn-l{z) (2.7) 
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P;;_i(z) = [S„(z)+v'(z)]P„_i(z)-A,_i(z)P„(z), (2i 



where 



An{z) := ^ f ""-^^^^^^^ P^{y)w{y)dy (2.9) 

l^n J a Z y 

Bn{z) ■■= j^ I "'^'^ y'^^^ Pn{y)Pn-i{y)w{y)dy. (2.10) 

A direct computation produces two fundamental supplementary (compatibility) conditions valid 
for all ^ G C U {00} and these are stated in the next lemma. 
Lemma 2 The functions An{z) and Bn{z) satisfy the conditions: 

Bn+l{z) + Bn{z) = {z- an)An{z) - v'(z) (Si) 

l + {z- an)[Bn+i{z) - Bn{z)] = /3„+iyl„+i - /?„A„_i(z) (5*2). 

It turns out that there is an equation which gives better insight into the coefficients a„ and 
Pn, if (Si) and (82) are suitably combined to produce a "sum rule" on An{z). We state this 
in the next lemma. The sum rule, we shall see later, gives important information about the 
logarithmic derivative of the Hankel determinant. 
Lemma 3 The functions An{z), Bn{z), and the sum 

n— 1 

satisfy the conditions: 

Bliz) + y'{z)Bn{z) + Y, M^) = Pn An{z) An-l{z). (S!,) 

j=0 
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Remark 5 If the v'(z) is a rational function in z then we observe that the divided difference 

z-y 

is also a rational function in z and y, from which we infer that (2.9) and (2.10) are rational 
functions of z. It is now clear that the compatibility conditions in Lemmas 2 and 3 would give 
further insights into the recurrence coefficients and certain auxiliary quantities that appear in 
the theory. 

2.2 Painleve V Continuous cr— Form and the Deformed Laguerre Weight 

In this subsection, we focus on the deformed Laguerre scenario, for which the main challenge is 
to characterize the numerator of (1.30). This is given by the following key result: 

Theorem 1: The Hankel determinant of the deformed Laguerre weight w{x) in (1.32) admits 
the following representation: 

Dn{t, A) = t"^ exp ( f ^^^M^^dx\ (2.11) 



where Hn{t) satisfies the Painleve V continuous Jimbo-Miwa-Okamoto c-form: 

{tH'lf = [tH'^, - Hn + K{2n + a + \) + n\] ^ - i{tH'^ - H^ + <5„) [{H'^f + A<] (2.12) 

with Sn '■= n{n + q + A). 

The remainder of this subsection is devoted to the proof of Theorem 1. 

2.2.1 Compatibihty Conditions, Recurrence Coefficients and Discrete Equations 

We start by noting that since 
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for the problem at hand is a rational function, the divided difference 

z-y 

will also be a rational function of z and y. Consequently, the functions An{z) and Bn{z) are 
rational in z. Therefore {Si) and (S'2) will give insight into various n and t dependent auxiliary 
quantities, which we shall see later. 

For the purpose of applying the ladder operator method, we carry out the preliminary 
computations, 



v(z) := — log w{z) = —alogz — X\og{z + t) + z 

V (z) = + 1 

^ ' z z + t 

v'(z)-v'(y) _ «^ A 



z-y zy {z + t){y + t) 

and we see that using the definition of our An{z), Bn{z) and applying integration by parts, 

An{z) = ^_l3M + Ml (2.13) 

„ / N n + rr,(t) r„(t) 

Bn{z) = ^^ + J^^ (2-14) 

Rnit) := ^ r ^-^^^w{y,t)dy (2.15) 

hn Jo y + t 

/in-i Jo y + t 

For the purpose of using (^i) and (S'2), in particular (6*2), we first state the following results 
obtained by substituting An{z) and Bn{z) given by (2.13)-(2.16): 

n-l 



Bl{z) + y'{z)Bn{z) + ^Aj{z) = z '^[{n + rnf + a{n + Vn)] 

j=0 

n~l 

+ z-^{n -Yj^:i + '■"[^ -a-t-2{n + r„)]/t + (n - A)/t} 



i=o 
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n-1 

j=0 
-2u2 



+ t)"^| Y^ Rj + rn[t + a - A + 2(n + r„)]/t - nA/t} 

j=0 



Now from (Si) we find, 



(2n + 1 + r„+i + r„) = q - a„(l - i?„) (2.f7) 

r„+i+rn = A-i?„(t + a„) (2.f8) 



and from (S'2) we find, 

(n + r„)2 + a(n + r„) = /?„(1 - i?„)(l - i?„_i) (2.19) 

n-J2Rj + Y[X-a-t-2{n + r„)] + "^ ^ ^^ = T ^^^ " ^«-i)^« + (^ " ^n-i)^n] (2-20) 
i=o 

5^ i?,- + y [i + a - A + 2(n + r„)] - !L = _/^ [(1 _ i?„)i?„_i + (l _ i?„„i)i?„] (2.21) 
i=o 

r2-Ar„ = /3„i?„i?„_i.(2.22) 

Remark 6 Observe that (2.20) and (2.21) are equivalent. We shall see later that (2.20), when 
combined with certain relations, performs the sum 

71-1 

automatically in closed form. 

This sum will provide an important link between the logarithmic derivative of the Hankel 
determinant with respect to t, Pn, and r„, which is an essential step in establishing the Painleve 
equation. 

Note that although the difference relations (2.17)-(2.19) and (2.22) look rather complicated, 
these can be manipulated to give us insight into the recurrence coefficients a„ and /?„. 

Now the sum of (2.17) and (2.18) gives us a simple expression for the recurrence coefficient 
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a„ in terms of i?„: 

a„ = 2n + l + a + A- ti?„. (2.23) 

From (2.19) and (2.22) have 

n{n + a) + Vnia + A + 2n) = /3„(1 - Rn - Rn-i) 

or equivalently 

/3n{Rn + ^n-i) = /3„ - n(n + a) - r„(a + A + 2n). (2.24) 

Now substituting (2.22) and (2.24) into either (2.20) or (2.21) to ehminate i?„ and Rn-i leaves 
us the fohowing very simple form for Yll=o ^i' which will play a crucial role later, 

n-l 

t^Rj =n{n + a + X)- (3n-trn. (2.25) 

j=0 

But in view of (2.23), we have 



n-l 


n-l 


l^R, - 


= n{n + a + A) — N| Oj 


i=o 


j=0 


= 


= n(n + a + A) + pi(n). 



(2.26) 
Comparing (2.25) with (2.26) gives 

Pi(") = -Pn -trn- 

Note that Pi(ra) also depends on t, although this is not always displayed. 

We are now in a position to find an expression for /3„ in terms of r„ and Rn- This is found 
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by eliminating Rn-i from (2.24) and (2.22) resulting in 



/3« 



1 



1 — Rn 



r„(2n + a + A) + 



Rn 



+ n{n + q) 



We summarize the above in the following lemma: 

Lemma 4 The recurrence coefficients an and /?„ are expressed in terms of the auxiliary 
quantities rn and Rn as: 



Furthermore, 



an = 2n + l + Q + A — ti?„ 
1 



/3n 



1 — Rn 



rn{2n + a + X) + 



• n ^" n 



Ar„ 



Rn 



+ n{n + a) 



n— 1 



n{n + a + X) - Pn- tr-„, 

Pn trn- 



(2.27) 
(2.28) 



(2.29) 
(2.30) 



2.2.2 t Evolution and Painleve V: Continuous <t Form 

In this next stage of the proof, we keep n fixed and vary t. The differential relations generated 
here when combined with the difference relations obtained previously will give us the desired 
Painleve equation. 

A straightforward computation shows that 



d^ 
dt 



log hn = Rn- 



(2.31) 



But, from (2.6), it follows that 



d^ 
dt 



l^n{Rn — Rn-l 
PnR'i 



n ~ ■^'n 



Rn 



(2.32) 
(2.33) 
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where the last equahty follows from (2.22). 
Differentiating 



0= / x°(x + i)^e-^Pn(x)P„-i(x)(ix 
/o 



with respect to t produces 



/"OO /'OO 

= A/ {x + tf-^e-'''Pn{x)Pn-i{x)dx+ x"(x + t)V^ 

JO Jo 



dt 



Pi(n)x" ^ + ... 



[^ Pn-l{x)Pn{x) . ., ^, d 

^/ ittt;^ ^t'(a;)ax + /i„_i — Pi(n), 

Jo 



x + t 



dt' 



Pn-l{x)dx 



finally resulting in 



dt 



PiW 



(2.34) 



Upon noting (2.4), this implies 



dan 

~dt 



^n+l 1"n- 



(2.35) 



Now differentiating (2.30) with respect to t and noting (2.34), we find 



dpi{n) _ df3n d 



dt 



dt dt 



(tr„ 



dvr. 






The above result combined with (2.33) gives 



dt dt Rn 



(2.36) 



Because (2.28) expresses /3„ as a quadratic in r„, we see that r„ satisfies a Riccatti equation; 



drn ri - Xr„ 



dt 



Rn 



Rn 
1 — Rn 



r„(2n + a + A) + 



"n ^''"n 



Rn 



+ nin + a] 



(2.37) 



31 



In fact there is another Riccati equation satisfied by Rn, which can be found as follows. Elimi- 
nating rn+i from (2.18) and (2.35), and upon referring to (2.28), we find that 

2rn = t—^ + A - Rnit + 2n + a + A - tRn). (2.38) 



Let us now eliminate r„(t) from (2.37) and (2.38), which results in a second order ordinary 
differential equation (o.d.e.) satisfied by Rn{t), where n, a, and A appear as parameters. A 
further linear fractional change of variable 

Rn{t) = l-- — oi y = l 



l-y{t) " 1-Rn{ty 

shows that y{t) satisfies a Painleve V: 

.."- 3y-l ^^/^2 y' I (y-l)' (»^y A^\ , {2n + l + a + X)y y{y + 1) ^ ^^.sg) 



^ 2y(y-l)^^^ t t^ \2'' 2yJ t 2{y - 1 

We note this is 

/a2 A2 , 1 

^IT' "y 2n + l + a + A, -- 

For the continuous (T-form of this Py , note that 

Hn ■■= t—logDn 
at 

, n— 1 

= t-Y.^ogh, 

i=o 



n-l 

'^3 



j=0 

= n(n + a + A) - /3„ - tr„ (2.40) 

= n(n + a + A) + Pi(n), (2.41) 

where the last two equations follow from (2.29) and (2.30) of Lemma 4. 
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From (2.34), (2.40), and (2.41), we obtain expressions for /3„ and r„ in terms of Hn and H'^, 

Pn = n{n + a + X)+tH'^-Hn (2.42) 

r„ = -K, (2.43) 



where ' denotes ^. 

All we need to do now is to eliminate Rn to obtain a functional equation satisfied by Hn, 



K and HI 



fiHn,H'n,H';,n,t) = 0. 



For this purpose, we examine two quadratic equations satisfied by Rn, one of which is simply a 
rearrangement of (2.28) and reads 

''"~^''" + (3nRn = (3n- r„(2n + a + A) - n(n + a). (2.44) 

-n-n 

The other follows from a derivative of (2.42) with respect to t and (2.33), 

/3ni?n-^%^ = tH'l (2.45) 

Solving for i?„ and l/Rn from the linear system (2.44) and (2.45), we find 

2l3nRn = (3n - rn{2n + a + X) - n{n + a) + tH'n 
'^C'^^ir^) " ^n-rn{2n + a + X)-n{n + a)-tH';, 

which we rewrite as follows 

= ^^ tH'; + it + 2n + a + X)H'n-Hn + nX 

tH'n-Hn + n{n + a + X) ^ ' 



2^ _ -tHJi + {t + 2n + a + X)H'n -Hn + nX 
Rn ~ [H'^Y + XH'^ • 

The product (2.46) and (2.47) gives us the desired continuous (j-form (2.12) 



(2.47) 
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It is finally worth noting that with 



Dn{t,X)=:&Dn, 



(2.48) 



then after a little computation we find that Dn satisfies the Toda moelcule equation [55] 






(2.49) 



2.3 Discrete cr— Form and the Deformed Laguerre Weight 

As an alternative to the continuous PV o"-form, the following theorem establishes a discrete 
cj-form satisfied by the logarithmic derivative of the Hankel determinant taken with respect t; 
which is a non-linear difference equation in n. 

Theorem 2: The Hankel determinant of the deformed Laguerre weight w{x) in (1.32) admits 
the representation (2.11), with Hn satisfying the discrete (T-form 



'n{n + a)t + {6'^Hn + t)[Hn - 5n] 



n2 



6'^Hn + 2n + a + \ + t 



A 



n{n + a)t + {5'^Hn + t)[Hn - 5„ 



n{n + a)t + {b'^Hn - t)[Hn - S, 



5'^Hn + 2n + a + \ + t 



S'^Hn + 2n + a + X + t 

Hn+l — Hn){Hn — -f^n-l) (2.50) 



where 



5 Hn '■— H, 



n+l 



Hn-l 



denotes a second order difference in Hn- 
The initial conditions are 



Hi{t) 



<" log D,{t,\), H2{t) = fAogD2{t,\) 



dt 



with 



Di{t, A) = fioit) , D2{t, A) = ^ioit)^i2it) - /i?(t), 
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and the moments are defined in (1.35). 

The proof foUows similar methods as in the continuous case in the previous subsection; 
namely, we express the auxiliary quantities Rn, r„ and the recurrence coefficient /?n in terms of 
Hn and Hn±i, and substitute these into (2.22). To begin, note that since 



Hn — 1 2_^ R 






we find that 



Hn+l — Hn — tRn 
Hn-l — Hn+l = t{Rn + Rn-l)- (2-51) 

Multiplying by /3„ and using (2.24), we obtain the following linear equation in /3„ and r„, 

(t + 5'^Hn)(3n - (2n + a + X)t r„ = n{n + a) t. (2.52) 

There is a further linear equation in /3n and r„, obtained by rearranging (2.25), which is 

/3n + trn = n{n + a + \)-Hn. (2.53) 

Solving for /3„ and r„ from (2.52) and (2.53) leaves, 

ft. = „(„, + . + A) - g.. + "'" + "" ^J^";- '''^" - ""; + " + ">' (2.54) 

n{n + a)t + {b'^Hn + t)\Hn -njn^a^ A)] 
*''" " 52^„ + 2n + a + A + t ' ^'^^^ 

The discrete ir-form is found by substituting (2.54), (2.55) and 

i Rn = Hn+l — Hn, 
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into (2.22), i.e., 



r^- Ar„ = PnRnRn-l 



2.4 Painleve IV Continuous cr— Form and the Deformed Jacobi Weight 

We now consider the Hankel determinant generated by the deformed Jacobi weight, 

x"i(l-x)"2-^(x + t)^, xG(0,l), (2.56) 

discussed in Section 1.3. In contrast to the deformed Laguerre weight, existing characterizations 
for this case are available. Specifically, such deformation was investigated by Magnus [51], 
where an auxiliary variable similar to our i?„ in the last section was found to satisfy a particular 
Painleve VI. The continuous cj-form associated with this Pyi was derived recently in [27]. Other 
related work dealing with this weight can be found in [56]. 

To state the results of [51] and [27] in our context, we must first introduce some additional 
notation. Let {Pra{x)} satisfy the orthogonality relations: 

/ P^{x)Pn{x)x''^{l-xr'-^{x + tfdy = hra{t)5^,n (2.57) 

Jo 
and 



hm(t) Jo 



Then 



yit) = l (l + ^)«-W 



2m + «! + 02 + 1 
satisfies the following Pyi'- 



y = o - + — 7 + ^—7 (y 
2 Vy y-i y + tj 

I i l_\ / 

t l + t y + t)^ 
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where 



1, ,2 a? (a2-A)2 1-A2 

i^i = -(2m + ai + 02 + 1) , 1^2 = -Y, 1^3 = ^ , ^^4 = — ^ — 



Furthermore, let 



H^{t) = -t{l + t)^^^^^^ + c,t + C2, (2.60) 



where 



ci := —m(m + ai+a2) (.^-olj 

C2 := -[2m(m + ai + 02) + (a2 — A)(qi + a2 — A) — A(ai - Q2 + A)]. (2.62) 

Then Hm satisfies the following cr-form of Pyi'- 

= {H'„,-rl){H'^-Tl){H'„,-rl){H'^-rl), (2.63) 



where 

ai + a2 — A a2 — A + ai 2m + ai + 02 — A 2m + ai + 02 + A 
n = ;; , ^2 = , rs = , r4 = , 



and «! and a2 are defined by (1.45). 

As for the single-user case, the Hankel determinant for the multi-user situation can also be 
expressed as an integral of H„i: 

T^ t^\ r^ r ^4.n\ f f Cl S + C2 - Hmis) , 

Dm{t) = -D„,ai,a2-AKac] * exp / ' ' ds 

\Joo •S(l + S) 
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Note here that we have tacitly assumed that t = Pi/{P — Pi) is strictly positive. However, 
we expect the equations (2.59) and (2.63) to be formally valid for all t € M. 

3 Characterization via the Coulomb Fluid Method 

In this section, we present an alternative characterization based on the Coulomb fluid method. 
As we will see, the key benefit of this approach is that it leads to simpler expressions than the 
exact results obtained via the ladder operator approach. Moreover, whilst this method is based 
on large-n considerations, we will show numerically that the approximations are very accurate 
for very small dimensions also. In fact, in the following section, we will demonstrate that the 
Coulomb fluid approach actually captures the exact distribution of the mutual information to 
leading order in n. 

Whilst the Coulomb fluid has been applied extensively in the context of statistical mechan- 
ics, it is relatively unfamiliar amongst the wireless communications and information theory 
communities. As such, in the following we will first present some basic background material, 
based mainly on [28], before deriving new results for both the single-user and multi-user MIMO 
systems of interest. 

3.1 Preliminaries of the Coulomb Fluid Method 

Consider a function of the form 

Zn{\) 



^n(O) 



exp[-(F„(A)-F„(0))] (3.64) 



where 



Zn{X) ■■= exp 

'(L,U)" 



-$(xi,...,x„)-A^/(5 



(3.65) 
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with 



$(xi,...,x„) := -2 ^ log\xj - Xk\ + n'^\/{xj) . (3.66) 

l<j<k<n j=l 

This expression embraces the moment generating function representations for both the single- 
user MIMO capacity (1.38)-(1.39) and multi-user MIMO capacity (1.57)-(1.59), with appropri- 
ate selection of the functions /(•) and v(-), and integration limits L and U. 

The key motivation for this representation is that it admits a simple intuitive interpretation 
in terms of statistical physics, as observed in the seminal papers by Dyson [29]. In particular, 
interpreting the eigenvalues xi,...,x„ as the positions of n identically charged particles, the 
function ^(xi, . . . , Xn) is recognized as the total energy of the repelling charged particles, con- 
fined by a common external potential n\/{x). The function f{x) acts as a perturbation to the 
system, effectively modifying the external potential. For sufficiently large n, we can approximate 
the particles as a continuous fluid with a certain (limiting) density, a{x), and assume that it is 
supported on a single interval (a, 6). This density will correspond to the equilibrium density of 
the fluid, obtained via the constrained minimization 

minF[o-] subject to / a{x)dx = I (3.67) 

Ja 

with 

/b ph ph 

a{x) [n y{x) + Xnf{x)) dx — n / / a{x) log \x — y\a{y)dxdy . (3.68) 



As a consequence of the Frostman Lemma [57, pg. 65], the equilibrium density satisfles the 
integral equation 

A f^ 

v(x) H — f(x) — 2 / In Ix — y\a(y)dy = A, x G [a, b] 
n Ja 

where A is the Lagrange multiplier which fixes the constraint that the equilibrium density has 
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total charge 1. See [57] for a detailed discussion. The above integral equation with logarithmic 
kernel is converted into a singular integral equation by taking a derivative with respect to x for 
X G {a,b), 



Ja x-y n 



where V denotes Cauchy principal value. 

If v(a;) is convex in a set of positive measure, the solution to this problem can be found [28], 
with the optimal cr(-) taking the form 

a{x) = (j(x, A) = cjo(x) + ^i^, (3.69) 

where 

^{h-x){x-a) f" y'{x)-y'{y) 

MV = ^-^ / s //. N. ^^dy (3.70 

2^ J a {x-y)^{b-y){y-a) 

denotes the limiting density of the original system (i.e., in the absence of any perturbation), and 
,(x. A) = A^(x) = —r^=A^==V t ^^^^^^^f'{y)dy (3.71) 



2T:'^yf{h-x){x-a) Ja V-x 

represents the deformation of this density caused by the external perturbation. The solution 
theory of singular integral equations can be found in the monographs [58], [59], and [60]. See 
also [61] for numerous examples on the application singular integral equations to problems in 
elasticity, and [62] for the Lp version of the theory. 

The boundary parameters a and h are chosen to satisfy the supplementary conditions 

'''^^^ dx = (3.72) 



a \/{h- x){x - a) 



40 



and 



xv'(x) 
27r J a ^J{h- x){x-a) 



dx = 1 



(3.73) 



With these results, for sufficiently large n, the ratio (3.64) is then approximated by 



Z(0) 



exp 



-A- 



qComm.{rp\ 
2^1 (i^ _ ^^Comm.(^) 



(3.74) 



where 



^Comm.^^) = / f{x)~0{x)dx, s}°'^'^-{T) = n [ f{x)ao{x)dx . (3.75) 

J a J a 

Remark 7: With the above results, the moment generating functions (1.37) and (1.56) 
become 



A^(A) w exp 



-A' 



cComm./'7^\ _ 



(3.76) 



which corresponds to a Gaussian distribution with mean and variance given by 



-Comm. 



/icouiomb = -S^°"'"'-iT) - nlogT, a2„,,„^b = -5i^°'^'^-(T) . 



-Comm. 



(3.77) 



Therefore, the outage probability can be obtained via 



-fout(C'( 



out J 



1 + erf 



Cout — ^ 

V2^ 



(3.78) 



Remark 8: Based on the Coulomb fluid method, the error exponent (1.18) is approximated 
as follows: 



f tpComm. (rp\ ] 

Er{R) = maxo<p<i <^ p^ ' ^ ^ ' - p(5p'^'^-(r) + nlog(T)) - pR 



(3.79) 
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The key challenge is to evaluate the quantities 5f'°'^'^'(T) and S^'^'^'i^^ for the single- 
user and multi-user MIMO scenarios. These problems are addressed in the following subsections. 
As we will see, in both cases we will need to solve numerous integrals which are quite complicated 
and are not readily available. Thus, to aid the reader, we have succinctly compiled the solutions 
to these integrals in the Appendix, along with some detailed derivations. 

3.2 Coulomb Fluid and the Deformed Laguerre Weight 

In this case, we have the particularizations 

/(x) = -log(r + a;), v(x) =x-/31ogx, L = 0, U = oo . (3.80) 

First consider the constants, a and h. These are determined by the equations. 



la ^/{b-x){x- a) 2vr J^ Y^(6-x)(x-a) 

With the integral identities (6.16)-(6.18), we obtain 



^ = l + f, V^b = f3, (3.81) 



which leads to 



a = 2 + /3-2Vl + /3, 6 = 2 + /3 + 2Vl + /?. (3.82) 

Now consider the limiting density, o"o(x). In this case, with (3.70), (3.81), and the integral 
identity (6.16), it can be easily verified that 

. s 1 \/(b — x)(x — a) 
aQ(x) = —-^ ^-^ ^, a<x<b (3.83 

ZTT X 
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which is the celebrated Marcenko-Pastur law [63,64]. Substituting this distribution along with 
f{x) = — log(T + x) into (3.75), and integrating using the identities (6.1), (6.3), and (6.5) gives 



52c°--(r) = -^ 



[a + b) log ' ' 



— va61og 



2 y 2 

^+V(T + a)(T + b))^-T^^ 

iV^+Vby 



(3.84) 



We note that an equivalent expression can also be obtained by changing variables x — t- 
(/3 + l)x and invoking an integral result from [65]. The derivation of our result here, based on 
applying the Schwinger paranietrization of the log function (6.15), has the advantage of being 
much more direct, and moreover it can be used to derive other integral expressions encountered 
with the Coulomb fluid approach, as shown in the Appendix. Such parametrization is ubiquitous 
in the analytical computation of integrals arising in quantum field theory; see for example [66]. 

For q{x), substituting f'{x) = —l/{T + x) into (3.71) and using the integral identity (6.10), 
we calculate 

«^) - , ,, ' . , (l - ^'^":'f"^' ) ■ (3.85) 

2^T^y{b-x){x- a) \ x + T J 

Substituting this into (3.75), and applying the integral identities (6.1) and (6.2) we find. 



■•Comm. 



5l-°'°'^-(T) = -21og 



+ 



2 \T + bJ 2 \T + a 



(3.86) 



Figure 1 plots the Gaussian approximation to the distribution of the normalized mutual 
information (per antenna) of single-user MIMO systems, based on combining (3.86), (3.84), and 
(3.77), and compares with the true distribution generated via numerical simulations. The results 
are shown for a relatively low SNR, P = 5 dB, and for various antenna configurations. In all 
cases the Gaussian approximation is very accurate, even for n as low as 2. The situation changes, 
however, when P is increased, as shown clearly in Figure 2. In particular, as P increases, it is 
evident that the distribution starts to deviate from Gaussian, and that this deviation is most 
significant for small n. This behavior is in line with the results of [67], which shows that the tail 
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+ Simulation 
Gaussian Approx (Coulomb) 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 l.S 

Normalized Mutual Information (nats/s/Hz/antenna) 




(a) Tlr = Tit = ' 



1.5 2 2.5 

Normalized Mutual Information (nat,s/s/Hz/antenna) 



(b) Ur = 2nt {nt = n) 



Figure 1: PDF of normalized mutual information /(x, y)/n for the single-user MIMO scenario 
(deformed Laguerre case). Results shown for SNR, P = 5 dB, and different antenna configura- 
tions. In all cases shown, the Coulomb fluid approximation is very accurate. 

of the distribution behaves like an exponential random variable, rather than a Gaussian, when 
P — 7- oo (and n fixed). We investigate this phenomenon in more detail in Section 4 (see also 
[36]). 

Figure 3 compares the Coulomb fluid approximation for the error exponent, based on com- 
bining (3.86), (3.84) and (3.79), with the true error exponent computed via numerical simulation 
of (1.21) and (1.18). Interestingly, we see that in all cases, including both low and high P, the 
Coulomb fluid approximation is extremely accurate. 

3.3 Coulomb Fluid and the Deformed Jacobi Weight 

In this case, we have the particularizations: 



/(x) = -log(^^), v(x) = -(/.! log x-v;2log(l-x), L = 0, [7 = 1. (3.87) 



First consider the limiting density ao{x). This can be obtained by direct evaluation of (3.70). 
Alternatively, we may apply a simple transformation of the limiting density />(•) corresponding 
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Gaussian Approx. (Coulomb) 
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Figure 2: PDF of normalized mutual information /(x, y)/n for the single-user MIMO scenario 
(deformed Laguerre case). Results shown for different SNR values. The Coulomb fluid approx- 
imation is very accurate when the SNR P is low, however it becomes less accurate (the mutual 
information distribution deviates from Gaussian) as P increases. 

to the classical Jacobi weight (1.48), given by [35] 



piy) 

Bn 
An^ A 
Bn ^ B 



vr 1 — y 



(3.88) 



1 



(2n + ai + a2 + 2)2 
1 



02 — «! — 'i\/n{n + Qi)(n + a2){n + ai + 02) (3.89) 



, .„ al- a\ + Ax/nin + ai)(n + a2)(n + ai + 02) (3.90) 

(2n + Qi + a2 + 2y L 

^V2 -Vi- 4V(1 + c/?i)(l + ip2){l + ipi+ ip2) 



{ifl +ip2 + 2)2 

1 
(991 + (^2 + 2)2 



c/?i - V?? + 4^(1 + m){l + ^2){l + ^1+ V2) 



(3.91) 
(3.92) 



Here, the quantities A and B were obtained by recalling that ai = rupi and 02 = nip2, and 
taking the limit n — t- 00. To relate this density to o"o(x), we first note that 



v(x) 



logw:iac{x) 



n 



(3.93) 
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Figure 3: Error exponent for the single- user MIMO scenario (deformed Laguerre case). Results 
shown for SNR, P = 10 dB. The Coulomb fluid approximation to the error exponent is very 
accurate for both low and high SNR ratios. 
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where wjacix) is the classical Jacobi weight (1.47). Now, from the relation (1.49) we have 



y'(x) = JacVjL 

nW3ac{x) 

*Jac(l - 2X) 



nw ja.ci'i- - 2a::) 
which after substituting into (3.70) gives the desired result 



ao{x) = -p{l - 2x) 
n 



(3.94) 



[1 + ((^1 + ^2)/2] ^^^ ^)(^ «) ^ a<x<b (3.95) 

IT X{l — X) 



with 



1-B , l-A 

a:=^-— , 6:=^-—. (3.96) 



Substituting this distribution along with f{x) = — log[(T + x)/(l — x)] into (3.75), applying the 
partial fraction decomposition 

V{b-x){x-a) _ 1 (^ ah ^ (l-a)(l-b) \ ^^^^^^ 



x{l-x) A/(6-x)(x-a) V X a;-l 

and integrating using the identities (6.1), (6.3)-(6.8), we find 

52C°--(r) _ olpgr ^/^^^ + ^/J^^ ^ ,/^iog^ i-(V^-^(i-a)(i-6))2 



n[l + ((/?! + V92)/2] ^ V\/T+^+^/rT^y V(^/a6+7(TT^)(TT6))2-r2^ 

' (T + 1)2 - ( v^(r + a)(r + 6) - v/(l - a)(l - &) )'^ 



+ V(l-a)(l-6)log 



4(l-a)(l-6) 



We note that an alternative solution was also computed in [10], requiring the numerical evalua- 
tion of a certain fixed-point equation. 

For q{x), we substitute /'(x) = — 1/(T + x) — 1/(1 — x) into (3.71) and apply the integral 
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Figure 4: PDF of normalized mutual information /(x, y)/n for the multi-user MIMO scenario 
(deformed Jacobi case). Results shown for n^ = 4 and n^ = 3. The Coulomb fluid approximation 
is very accurate when the signal-to-interference ratio P/Pi is low, however it becomes less 
accurate (the mutual information distribution deviates from Gaussian) as P/Pi increases. 



identity (6.10) along with (6.11) to yield 



~e{T) 



27rY^(6 — x){x — a) 



V(l-a)(l-6) ^{T + a){T + h) 



1 — X 



x + T 



(3.98) 



Substituting this into (3.75), and applying the integral identities (6.2), (6.5), (6.8), (6.9), we 
obtain 



5r'^'"-(r) = iog 



4V(r + a)(r + 6)(l-a)(l-6) 



(T + 1)2 - (^(1 - a)(l - 6) - V(T + a){T + h) 



(3.99) 



Figure 4 plots the Gaussian approximation to the distribution of the normalized mutual 
information (per antenna) of multi-user MIMO systems, based on combining (3.99), (3.98), and 
(3.77), and compares with the true distribution generated via numerical simulations. As for the 
single-user MIMO scenario, the Gaussian approximation is accurate regardless of the number 
of interferers iC, especially when the signal-to-interference ratio P/Pj is not high. As P/ Pi is 
increased, once again the distribution starts to deviate from Gaussian; however, this deviation 
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Figure 5: Error exponent for the multi-user MIMO scenario (deformed Jacobi case). Results 
shown for n^ = 4 and n^ = 3. Similar to the single-user MIMO case, the Coulomb fluid 
approximation to the error exponent of multi-user MIMO channels is very accurate for both low 
and high signal-to-interference ratios. 

is seemingly less significant than that seen previously for the single- user case. 

Figure 5 compares the Coulomb fluid approximation for the error exponent, based on combin- 
ing (3.99), (3.98) and (3.79), with the true error exponent computed via numerical simulation of 
(1.21) and (1.18). Again, we see that in all cases, the Coulomb fluid approximation is extremely 
accurate. 

4 Beyond the Coulomb Fluid Approximation 

In this section, we take a closer look at the Painleve and Coulomb fluid representations of the 
mutual information distribution. We will focus on the Laguerre scenario (i.e., the single- user 
MIMO case), although the analysis can be extended to the Jacobi scenario also. Our main ob- 
jective is to establish relationships between the exact characterization of the mutual information 
distribution via the Painleve equation, and the simpler Coulomb fluid approximation. As a key 
result, we will show that for both the mean and the variance, the Coulomb fluid gives an exact 
representation to leading order in n. We also employ the Painleve equation to establish the first- 
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order correction terms to the mean, variance, and third cumulant, which helps to characterize 
the deviation of the mutual information distribution from Gaussian. 

4.1 Initial Study using Power Series Expansion 

To get an intuitive feel for the problem, we start by looking for a 1/t expansion in the Painleve 
equation as t — t- oo. This is essentially the case where the SNR, P, is small. For simplicity, we 
consider the case m = n; although the analysis can be extended to the case m^^ n. 

4.1.1 Evaluating Cumulants from the Painleve 

To proceed, let 

iin{t)=-n\^Gn(t) (4.100) 



such that 



M{X) = e^p(f ^^^dx) . (4.101) 



Note also that 



tj^logM{X) = Gn{t). (4.102) 

From the Painleve equation (2.12), Gn{t) satisfies 

{tG'^^f = {G'^it + 2n + X)-Gnf- 4(tG'„ -Gn + n^){G'^ + AG'„) (4.103) 

where the derivatives are with respect to t. Suppose that Gn{t) has a formal power series 
expansion in 1/t: 



oo 



Gn{t) _ V^ ^ 

rfiX ~ ^¥ 
k=l 
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(4.104) 



where the coefficient bk depends on n and A. As a first step, we substitute the power series into 
(4.103) and find the first few ft^'s as foUows: 

61 = -1 

62 = 2n - A 

63 = _(i + 5n2_6nA + A2) 

64 = lOn + Un^ - (5 + 29r?)X + 12nX^ - X^ 

65 = _(8 + 70n2 + 42n*^-(80n + 130n^)A + (15 + 95n2)A2-20nA^ + A^) 

66 = 168n + 420n^ + 132n^ - (84 + 80671^ + 562n'')A + (350n + 624n^)A2 

-(35 + 235n2)A^ + 30nA^ - A^ 

67 = -{180 + 2121n2 + 231071"^ + 429n^ - (2366n + 6510n3 + 2380n^) A 

+ (469 + 4795n2 + 3682n*^)A2 - (1120n + 2128n^)A^ + (70 + 490n2)A^ - 42nA^ + A^}. 

From these examples, it is clear that bk takes the form 



fc-i 



6. = E^1^ 



,fc-^-i 



Bkin'^'-' + Ckm'-'~' + 0(n 



k-e~3 I n/^A;-£-5> 



(4.105) 



e=o 



where the coefficients B^/ and Ck/ are independent of n and A, and are computed via 



Bk,£ 



d 



•k-i-l 



d' 



{k - I - l)U\ dn''-^-^ dX^ 



bk 



A=0,n=0 



A; = 1, 2, ... , 00, 



0,...,A;- 1 (4.106) 



and 



Cfcj 



1 d^-^-3 d' 

{k - e - sy.ei dn''-^-^ dX^ 



bk 



A=0,n=0 



A; = 3, 4, ...,00, e = 0,...,k-3 (4.107) 



respectively. Note that Ck/ = for k = l,k = 2, oTi>k — 2. The B^./ coefficients represent 
the leading order terms in n, whereas the Ck/ coefficients represent the next (non-zero) lower 
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order correction terms; for example, 

Bifl = —1, -62,0 = 2, i?3,o = -5, i?2,i = -1, -83,1 = 6, i?3,2 = -1 

and 

C3,o = -l, C4,o = 10, C5,o = -70, C3,i=0, C4,i = -5, C5,i = 80. 

We aim to investigate the cumulants of the mutual information distribution, which in turn 
requires an expansion of the form 

Gnit) = Xgi{t) + X'g2it) + --- (4.108) 

Together with (4.104) and (4.105), we then get 

gk{t) = gk{n/P) = n"-^ f; (.B.^k-i + ^^ + O (^] ) P^ k = l,2,... (4.109) 

Plugging this into (4.101) and integrating, the logarithm of the moment generating function 
becomes 

logA^(A) = J]Kfc— (4.110) 

fc=i 



where k^ denotes the /cth cumulant given by 



,, [^ 9k{n/y) 

Kk = -k\ / dy 

Jo y 

e=k 



The mean is 
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the variance is 



., . ,,_,,,^^,._^^^,„^j_a,._,,^,s^,.,^,oa 



the third cumulant is 



K3 



n 



1 



IP" - 18P^ + 1UP° + • • • + ^ ilSP" - 350P*' + • • • 



6 



n^ 



+ 



w 



and so on. 

Importantly, this result demonstrates that as n grows large, the kth. cumulant scales as 



Kk = 0{n'-^), A; = 1,2,... 



(4.112) 



Thus, as expected, we see that the mean of the mutual information grows linearly with n, the 
variance converges to a constant, and all other cumulants disappear as n — )■ oo. This reaffirms 
that the distribution becomes Gaussian for asymptotically large n. 

Whilst in the analysis above we have substituted the numerical values of the constants 
Bij., explicit formulae can also be derived by directly using the Painleve differential equation. 
We demonstrate the procedure by considering the coefficients in the mean summation (4.112), 
namely {B^^q}. The same technique can be used to derive the formulae for the coefficients of 
the higher moments also. 

Plugging the series (4.104) into the differential equation (4.103), keeping only the lowest 
order terms in A and taking n large, the r.h.s. of (4.103) becomes 



{PnXf 



oo 



(4^2 + 45i,o) + (1251,0^2,0 + 8^2 + 8B2,o) P + ^ Xt+2P 



i=2 



(4.113) 



with 



j-i 



Xi+2 = -4S,+i,o - 8iB^fi + J^(2 + j)Bj+i,o {{i -j + 2)S,_^+i,o + 4(i - i)Bi_^-,o) • (4.114) 
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Now consider the l.h.s. of (4.103). We have 



(«i;f^^f,.(, + l)B..oP-^-V^„ (4,115) 



(nA)2 \~~f '^ 



j=i 



as n — 7- oo. 



Thus, dividing both the l.h.s. and r.h.s. by (PnA)^ and taking n large we get 



oo 



{ABIq + 45i,o) + (l2Si,o52,o + 8B(o + 8^2,0) P + 2_^ Xt+2P' = . (4.116) 

i=2 

We calculate the Bi^Q's recursively, since every coefficient of P must equate to zero. Trivially, 
considering the constant and linear terms in P, we get 

5i,o = -1, B2,o = 2 . (4.117) 

For higher order coefficients, we have the recurrence relation 

Bi+i^o = -2iB,fl + - ^(2 + i)B,+i,o {{i -3 + 2)5,-j+i,o + 4(i - i)S,„,, o) (4.118) 

i=i 

for i > 2. Based on this, the next few coefficients are evaluated as 

53,0 = -5, ^4,0 = 14, ^5,0 = -42. (4.119) 

It can be verified that this difference equation also admits the explicit solution 

and therefore the mean of the mutual information takes the explicit form (to leading order of n) 
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4.1.2 Comparison with the Coulomb Fluid 

Now, consider the corresponding quantity derived based on tlie Coulomb fluid method: 

Gn{l/P) = tj^logMiX) 

= X~gi{l/P) + X^g2{l/P) (4.122) 



where 



d 



~9i{l/P) = P^{S^'''^'^-{l/P) + nlogP) 

UVP) - P^'-?^^:^^ . (4.123) 

Here we have substituted the expression for A^(A) given in (3.76), and used the fact that t = n/P. 
From (3.84), we compute 

5i(l/P) = -n 1 + ^ (l - Vl + 4P) (4.124) 

which, after applying a Taylor expansion of \/l + AP around zero, gives 

oo 

~gi{l/P) = nY,Bk,oP'' . (4.125) 

k=l 

To leading order in n, this agrees precisely with gi{t) in (4.109), thereby establishing that the 
Coulomb fluid method gives the exact value of the mean mutual information for large n. 
From (3.86), we compute 

S^°'^"'-{l/P) = -P2 + 4p3 _ _^p4 ^ 52P5 - ^^p6 ^ ggQpT (4^26) 

giving 

g2{l/P) = -P^ + 6P^ - 29P'' + 130P^ - 562P^ + • • • (4.127) 
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This agrees precisely with (72 (i) in (4.109) to leading order in n, thereby confirming that the 
Coulomb fluid gives the correct asymptotic variance, in addition to the correct asymptotic mean. 
In summary, we have 

logX(A) = ^(^Co.,„.. + iQp»4p* + fp» + -)+o(J3)) 

\3 / 1 1 /I 

+ — - (2P3 _ i8p4 + ii4p5 + . . \ + CispS _ 35Qp6 ^..\^qI 
3! \n ^ ' n-^ ^ ' \n^ 



fc! V ^ 

Here, /xcouiomb and cr^f^^if^^^ represent the mean and variance respectively, calculated based on 
the Coulomb fluid method. All other terms represent correction terms, which essentially account 
for the deviation of the mutual information distribution from Gaussian for flnite values of n. 



4.2 Refined Analysis for All P 

We now present a more refined analysis, which does not require a power series representation 
of 1/t. This analysis is based on evaluating a non-perturbative summation of the perturbation 
series in P, obtained from the terms of Gn{t) which are linear in A, and to leading order in n. 
For this purpose, we substitute the series representation (4.108) into (4.103). 

4.2.1 Analysis of the Mean 

We start by considering the first cumulant (i.e., the mean). By comparing the coefficients of 
A on the l.h.s. and r.h.s., we find that the coefficient of A is identically equal to 0, and the 
coefficient of A^ satisfies: 

(51)' - An'g^ - 2[t + 2n]gig'^ + [t" + Ant]{g[f - t\g'if = 0. 
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Note that 



t = — =: nT. 
P 



We have introduced T so that the differential equation does not get too comphcated. After this 
change of variable in t to nT, without introducing further notation in place of (71, the differential 
equation becomes 

(<7i)2 - ^ng[ - 2(r + 2)g^g[ + [T^ + 4T]{g[f - ^"^ = 

n^ 

where now ' denotes the derivative with respect to T. Letting 

gi{nT)=nY(J) (4.129) 

we find that 

y2 _ 4y/ _ 2(T + 2)YY' + [T^ + 4Tly'2 ^-^ = 0. {Y) 

n^ 

Taking n — )■ 00, y is seen to satisfy 

y2 _ 4y/ _ 2(7^ + 2)yy' + [T^ + 4r]y'2 = o. (yo) 

Now consider gi{T), the term analogous to gi{nT) but derived based on the Coulomb fluid. 
From (4.124), 

UT) = nY^{T) (4.130) 

where 

Yo{T) = ^, ^ ^ . 4.131 

4 + T + VT(4 + T) ^ ^ 

This expression is found to satisfy Eq. (YO) identically, thus confirming that the Coulomb fluid 
approach gives the exact value for the mean of the mutual information to leading order in n for 
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all values of P. Note that this result is stronger than that derived in the previous section, since 
it appHes even for values of P (or values of 1/t) for which a formal convergent power series in 
1/t does not exist. 

To compute the 1/n correction to the previously obtained Y, we substitute 



Y{T) = YoiT) + \Yi{T) + 0^ ^ 



n^ 



n^ 



(4.132) 



into Eq. (Y), and then obtain Yi by setting the coefficient of 1/n^ equal to 0. This gives 

2^0^! - 2(2 + T)YiY^ - AYl - 2(2 + T)YqY{ + 2T(4 + T)YqYI - T^ Yq = 

which is a first order linear equation in Yi. However, the coefficient of Yl vanishes identically 
when we make use of Yq from (4.131). The solution of the algebraic equation reads: 



Yi{T) 



T(4 + r)5/2 



With these results, we can compute the asymptotic mean of the mutual information, includ- 
ing the first-order corrections, as 



Kl 



9i{n/z) 



dz 



/^Coulomb ~r /^Correction ~r v-' 



n 



n^ 



(4.133) 



where jUcouiomb is the mean value computed via the Coulomb fluid, which from (3.84) and (3.77) 
with /3 = is given by 



/^Coulomb — ^ 



2 log 



1 + Vl + 4P\ (1 - VI + 4P) 



4P 



and /Ucorrcction is the first order correction term given by 



/^Correction 



yi{i/^) 



dz 



1 

12 



1 + QP + 6P2 

(4P + 1)3/2 



(4.134) 



(4.135) 
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Note that if we expand this expression for ki around P = 0, the series matches precisely with 
(4.111) as expected. 

4.2.2 Analysis of the Variance 

Now we consider the second cumulant. To this end, after substituting (4.108) into (4.103) and 
setting the coefficient of A^ equal to 0, we get 

25^(1 + g[){gi - tg[) - 2n\l + 2g[)g'^ 

+{-gi + (2n + t)g[}{-g2 + g[ + (2n + t)(7^} - t^g'l g'^ = 0. (4.136) 

Once again, applying the change of variable t = nT, but without introducing new notation for 
gi and g2, we find 



-g[(l + '-^) {gi-Tg[) -2n(l + 2'-!^) g'2 



n ^ ^ 

1 



+{-51 + (T + 2)5'i}{-52 + g'xln + (2 + r)<7^} - ^T^ g'{g'i = (4.137) 



with ' denoting d/dT. 
Let 



52(nT) = Zo(r) + \zi{T) + o(\] . (4.138) 



n^ \n 



We first compute Zq{T), the leading order term in n. To this end, substituting (4.138) along 
with (4.129) and (4.132) into (4.137), and then keeping only the leading order terms in n (the 
terms which are linear in n), we obtain 

2yo'(l + Yi){Yo - T FoO - 2(1 + 2Y^)Z', 

+{-Yo + (T + 2)Y^}{-Zo + (T + 2)4 + Y^} = 0. (4.139) 

Interestingly, if we plug in the expression for Yq{T) given in (4.131), we find that the coefficient 
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of Zq is identically equal to 0. Thus, (4.139) reduces to a simple algebraic equation, whose 
solution is: 



2 2 V 4 + r 4 + r 

After making the substitution T = 1/P, it can be verified that this expression matches precisely 
with g2{T) in (4.123), derived based on the Coulomb fluid method. This result confirms that 
the Coulomb fluid approach gives the exact value for the variance of the mutual information 
to leading order in n for all values of P. Again, this result is stronger than that derived in the 
previous section, since it applies for values of P (or values of 1/t) for which a formal convergent 
power series in 1/t does not exist. 

Now consider the correction term, Zi{T), in (4.138). Again we substitute (4.138) along with 
(4.129) and (4.132) into (4.137). In this case, however, we extract only the terms of order 1/n, 
which gives a rather large first order equation in the unknown Zi(T). Fortunately, we find that 
by plugging in the previously determined equations for Yq (T) , Yi (T) , and Zq (T) , the coefficient 
of Z[{T) vanishes identically, and so we are left with a linear equation in Zi(T). This is easily 
solved and we find 

^ ,^, 8 + 16T + 20T2 + 5T^ + 6v^T(4 + T) + WT^/'^Vi + T + hT^I'^^J\ + T 

Zi{T) = y , ■ — N • (4.140) 

r(4 + r)4 f 2 + 4r + T2 + 2y^T{4 + T) + r3/2 V4TTj 

With these results, we can compute the asymptotic variance of the mutual information, 
including the first-order corrections, as 

r92in/zl, 
K2 = — 2 / dz 

Jo z 

= '^Coulomb + ~2 '^Correction + O ( "4 j (4.141) 

where cr^ouiomb ^^ ^^^ variance computed via the Coulomb fluid, which from (3.86) and (3.77) 
with /3 = is given by 



•^Coulomb — 2 log 



(4P + 1)1/4 + (4p + l)-l/4 
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(4.142) 



and <7^orrection ^^ ^^^ ^^^^ Order correction term given by 



'^Correction 



Zi{l/z] 



dz 



1 

12 



8P3(1 - 3P) 12P3 + 30p2 ^ iQp _^ I 



(4P + 1)3 



(4P + 1)5/2 



. (4.143) 



Note that if we expand this expression for K2 around P = we get a series which matches 
precisely with (4.111), as expected. 



4.2.3 Analysis of the Third Cumulant 

Now consider the third cumulant. After substituting (4.108) into (4.103), setting the coefficient 
of A'^ equal to 0, and then going through the same procedure as before (i.e., applying t = nT) 
we get 



4 1 + 2 



(51 



n n \ n J 



+ (-52 + fj- + (2 + T)g'^ + 2(-5i + (T + 2)g'^) (-33 + f + (2 + T)g'^ 



-4n^ 



^V+fl + 2^')^ 
n I \ n I n 



w 



-{g2+'^9ig3) =0 



(4.144) 



where the derivatives are taken with respect to T. Now substitute (4.138) along with (4.129) 
and (4.132) and 



53(T) = ^Ao(r) + -j3Ai(T) + oQ 



n 



n^ 



n" 



(4.145) 



where Xq{T) and Ai(T) are n independent. We find that the highest order term in n is in fact 
n independent, and using the previously determined equations for Yq{T), Yi(T), Zq(T), and 
Zi (T) , once again we find that the coefficient of X'q (T) vanishes identically. The resulting linear 
equation in Ao(T) has solution 



(4.146) 



, ._l ( 2\/r T \/T \ 

0^ >-2 I (4 + T)5/2 + (4 + T)2 " (4 + r)3/2 1 
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Proceeding in a similar manner, we find 

16 30 aVT 9T3/2 ^.5/2 

1*^ )- ~ 7^3/2(4 _^ 2^)11/2 ~ 7^1/2(4 _^ 7^)11/2 ~ (4 + r)iV2 + (4 + ^)11/2 ~ 2(4 + r)iV2 

8 10 lOT T^ 

+ r(4 + T)5 + (4 + T)5 (4 + r)5 ^ 2(4 + r)5 ' ^ ^ ^ 

With these results, we can compute the asymptotic third cumulant of the mutual information 



as 



.3 = -31 r '-^^^^^dz 

Jo z 



— /^S, Correction A H ^/^S.CorroctionB + O ~^ (4.148) 



where 



;^Xo(iA), 

/tSjCorrectionA — "" / "^ 



z 

^ 5^ + - (4 1491 

^TTIP l+4P^ (l + 4P)3/2 ^- ^ 



and 



/-^XiO/z) 

/^S, Correction B — " D / "2: 

Jo 2 

16P6 48P^ 380P^ 610P'' 323P^ 



(4P + 1)4 (4P+1)5 3(4P + 1)5 3(4P + 1)5 3(4P + 1)5 
80P2 lOP 1 42P5 93P'^ 

~ 3(4P + 1)5 ~ 3(4P + 1)5 ~ 6(4P + 1)5 + (4P + 1)9/2 + (4p + 1)9/2 
71P3 21P2 3P 1 

+ (4P + 1)9/2 + (4P+ 1)9/2 + (4P+ 1)9/2 + 6(4P + 1)9/2 ' ^ '^^^^ 

Note that if we expand this around P = we get a series which matches precisely with (4.111). 
The correction terms to the mean, variance, and third cumulant are plotted in Figure 6. In 
particular, the "Mean" curves represent 

n(£'[/(x,y)] - /^Correction) 
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with /Ucorrcction given by (4.135); the "Variance" curves represent 

n2(Var[/(x,y)] - a^orrcction) 

with cr^orrcction gi^en by (4.143); and the "K3" curves represent 

n(K3[/(x,y)] - Ks.CorrcctionA) 

with K3, Correction A given by (4.149). Here, the mean E[I{x,y)], variance Var[/(x, y)], and third 
cumulant K3[/(x, y)] of the mutual information were calculated using numerical integration 
procedures in Maple. From the figure, we can make some interesting observations. First, it is 
clearly evident that for low SNR P, all three correction terms converge very quickly to zero, 
confirming the near-Gaussian behavior of the distribution even for small n, as seen previously 
in Figure 2. However, as P increases, all three correction terms become much more significant. 
This, in turn, leads to a larger deviation from Gaussian, which again is in line with the numerics 
presented in Figure 2. It is also particularly interesting to note that the correction terms for the 
higher-order moments tend to deviate quicker than the lower order moments. To understand 
this phenomenon, it is useful to look closer at the correction terms as P grows large. This is the 
focus of the next subsection. 

4.3 Analysis at Large P 

For large P, the asymptotic mean (4.133), variance (4.141), and third cumulant (4.148) become 

^ / 1 



^ ~ n log P + 

n \lo 

2 log(2P) , P fS 
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Figure 6: Correction terms to the mean and variance of the mutual information for the single- 
user MIMO scenario (deformed Laguerre case). Results are shown for rir = rit = n. 
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From these expressions, we notice that in all three cases, the correction terms (i.e., the second) 
are increasing in P. Since these terms give corrections to the Coulomb fluid Gaussian approxi- 
mation, this implies that as P increases, the Coulomb fluid approximation loses accuracy, and 
the distribution deviates away from Gaussian, as our previous numerical results have indicated. 
We can also gain insights into the relative effect of P on each of the correction terms. Specifi- 
cally, for the mean, we see that for the correction term to dominate the leading term, it must be 
at least 0(n'^). The variance, on the other hand, must only be O(n^), whist the third cumulant 
is even more sensitive to P, and only requires P to be 0{n). This conflrms that the higher 
moments are more sensitive to the variations in P, and moreover, it allows insights into the 
Gaussianity of the mutual information distribution in terms of both n and P. 

We would like to mention, however, that some caution should be exercised in interpreting 
these results. In particular, since our analytical results which are asymptotic expansions in n are 
based on the Painleve V representation of the moment generating function, they are formally 
valid for large n but fixed T or P. Therefore, it is to be expected that when P increases to the 
point where the corrections are comparable with and/or overwhelm the main terms, the Coulomb 
fluid approximation will break down. Providing a rigorous investigation of the situation where 
both n and P increase simultaneously is the subject of on-going work. 

4.4 Asymptotic Recurrence CoefRcients 

In addition to deriving the asymptotic moments, having determined the large n expansion of 
GninT), we can also compute the recurrence coefficients, a„ and /3„ for a = and large n with 
t = nT. For this purpose, note that a„ is easily related to i?„ (see Eq. (2.20)) and i?„ is in 
turn expressed in terms of Hn (or equivalently G„) and its derivatives (see Eq. (2.39)), while 
Pn related to Hn and its derivative through Eq. (2.35). 

From our previous analysis, we have learned that for t = riT and n large, Gn{nT) has the 
expansion 

Gn{nT) = x(nYo{T) + ^^ + ..)+X'(zo{T) + ^^ + ... 



65 






From this series, after a straightforward if lengthy computation we find that 



an{nT) = (2-^)n + A(l 



(4 + T)2 
and 



2/ \ V 4 + T 

+ A- 



VI+T + X{V4 + T - Vt)] - + O ( \] , (4.152) 

J n yw^ J 



j3n{nT) = n^ -\ , n 



2T + (2 + r)(4 + T)(-l + VT/^/4TT) 



2(4 + T)2 

^^ 6 + 4r + A^(2r=^ + r3-r5/2^/4TT + 2T) i qM > .^^^g^ 



5 Concluding Remarks 

This paper has focused on the computation of Hankel determinants which arise in the information- 
theoretic study of MIMO communication systems. We considered two practical scenarios; 
namely, single- user systems, in which case the determinants of interest are generated from a 
certain deformed Laguerre weight, and multi-user systems, in which case they are generated 
from a deformed Jacobi weight. In both cases, we derived new exact characterizations of the 
Hankel determinants in terms of classical Painleve differential equations, as well as closed- form 
expressions which are formally valid for large matrix dimensions, but give accurate results for 
small dimensions also. 

We also demonstrated that, in contrast to most prior work dealing with MIMO information 
theory, our exact and approximate results can be used together to derive explicit formula for 
the mean, variance, and higher order cumulants of the mutual information, as well as their 
corresponding correction terms. This is important, since it allows one to gain insights into the 
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Gaussianity of the mutual information distribution in terms of key system parameters, such as 
the number of antennas and the signal-to-noise ratio. For example, by focusing on the single- 
user MIMO scenario, we showed that the mutual information distribution may deviate strongly 
from Gaussian when the signal-to-noise ratio is relatively high, and we explicitly captured this 
effect through the cumulant formulas. 



6 Appendix: Some Relevant Integral Identities 

For the Coulomb fluid derivations, we will require the following integrals: 



:^M^±^d. = 2.1ogf^^^ + ^^^l (6.1 

' l"g(^ + ^) dx = ^-^ loJ^— + ^—] (6 2 

a ^/{b-x){x-a){x + t) ^{t + a){t + b) ^\2Vt + ^ 2^/t + h ) 



; iog(x+t^,, = .Uj^>W¥±^p3?^\ (6.3 

a ^J {h - x){x - a) x Vab \ i\^ + Vby J 



^ xlog{x + t) ^ (VoTt - \/6Tt)2 
ax = TT 



a \/{b- x){x - a) 



^^a + b ^^^ ( {^{a + t){b+^t)+tf _-ah \ ^^^^ 



f 

J a 




log(a: 


+ t) 




r- 


x)(x - 
' lo 


-a)ix- 
g(l-x) 


1) 




Vib 
log 


- x){x — 
(1-x) 


a) 


t 


log(l 


x){x — a 
-x) 


)x 


L 
t 


Vib- 


x){x - 
log(l 


- a) {x — 
-x) 


1) 



log 



(v^I^+v^I^)^ 



V {t+l)2-{ ^/{t+a){t+b)-^(l-a)(l-b) Y 

ax = vr , 6.5 

V(l-a)(l-6) 

dx = 2.1og( ^^^^ + ^^^^ ) (6.6 



TT 



l-(V^-V(l-a)(l-&)) 



2 



dx = -^= log ^ ^/ ^ JK J_J_ (g.7 

Vab ^\ (^ + V6)2 ) ^ 

27r , / 1 1 \ 

dx = = log , + , (6.8 

v^(i-a)(i-5) V2\/T^^ 2^/r^6y ^ 

,' it+lf-iy/it+a)it+byy/(l-a){l-b)f 



(Vt+^+Vt+b)^ 

- (^'j^ ^ j^ -^^ -^ 1— (69 

la yj{b-x){x-a){x + t) y/{t + a){t + b) 
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Before proving these results, we state for reference the following identities: 

dx 1 ^ 2v/c\/ax2 + bx + c + bx + 2c „ /^ , „x 

log , OO, (6.12) 



a;\/aa;2 + 6x + c \/c x 

/dx 1 r 
= -^ log 2Vaby/{t + a)(t + 6) + 2a6x + (a + b)t (6.13) 
y (t + ax)(t + 6a;) va6 '- -' 

dx _ o + ^ 1 /(\/(a: + a)(x + 6) + y/abf - x'^\ 

x^y/{x + a){x + Vj ~ 2(a6)3/2 °^ I x ) 



V^(x + a)(x + 6) 



(6.14) 



abx 
log{A + B) = \ogA + j -^^ (6.15) 

^" " (6.16) 



a (x + t)v^(6-x)(x-a) v^(t + a)(t + 6) 



a Y^(6 — x)(x - a) 



vr (6.17) 



xdx a + 6 

' ~~f^^^=^^^^ = ^ — ^ — (6.18) 

a ^J{b-x){x-a) 2 

V I . "^^ = (6.19) 

Ja {x-y)^J{b-y){y-a) 

We start with some brief remarks concerning the identities (6.12)-(6.19), before focusing on 
the proofs of the main results (6.1)-(6.11). The integral (6.12) is [68, Eq. (2.266)], (6.13) is a 
minor variation of [68, Eq. (2.261)], and (6.14) follows from [68, Eq. (2.269.2)]. Note also that 
the r.h.s. of equation (6.15) is 
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which obviously equals the l.h.s. To give an indication how (6.16)-(6.19) may be proved, first 
consider the analytic function 

R{z) = y/{z-a){z-b) 

defined in the complex plane slit along [a,b]. Here we assume that 

< a < 6 < 1 
without loss of generality. The branch of R{z) is chosen in such a way that 

R{z) -^ z, as K z — ;• oo. 

Let F{x) be defined for x G M and extended to F{z), a meromorphic function with poles. Let 
A be a "dog bone" contour traversed clockwise above and below the segment [a,b], where the 
point of oo is contained in the interior of A. Keeping in mind that 

R±{x) = ± jy/{b - x){x - a), X £ {a,b), 

where R±{x) is defined to be the analytic continuation of R{z) to above and below the segment 

(a, b), we have that 

F{x)dx j f F{z)dz 



la \/{b — x){x - a) 2 y^ ^(^z - a){z - b) 
and the r.h.s. of the above equation can evaluated using residue calculus. The equation (6.16) 
follows immediately by computing the residue at t, while in computing (6.17) and (6.18) we 
should keep in mind the contributions from the residues at oo. To compute the principal value 
integral (6.19), we first define 

gyt) :- 



a {y + t)^{b-y){y-a) ^{t + a){t + b)' 
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The principal value integral is then evaluated as 

-olimb(-a; + ie) +g{-x-je)], X G (a, 6). 

An easy computation gives (6.19). 

We now come to the main integrals, (6.1)-(6.11). Start by considering (6.1)-(6.5). Of these, 
we will explicitly derive (6.2); the other integrals are evaluated in a similar way with the help 
of the properties (6.12)-(6.14). Using (6.15) along with (6.16), we obtain 



J a 



log(a: + t)dx 



{x + t)y^{b — x){x — a) 



V(t + a)(t + 6) ^^io A-1 l^^(t + a)(t + 6) ^ [t + a\){t + h\) j 

^^°^^ +.C'-i , ' -^ ' ]dx 



^/{t + a){t + b) Jo X y ^{-ax + t + a){-bx + t + b) ^^{1 + a){t + b) ^ 

+ vr lim / — — , dx, 



y^{t + a){t + b) ^^oJe X y ^{-ax + t + a){-bx + t + b) ^J {t + a){t + b) ^ 

where we have made the substitution x = 1 — A and have replaced L ... by linie^o ji\ •.., so that 
we may invoke (6.12). The integration is now completed as 



J a 



log(x + t) _ vr log t 



{x + t)^{b-x){x-a) ^/{t + a){t + b) 

+ E"o vc + »)« + 1) (" '"" ^e + ")(« + ") " '°^"'^" + '■)(' + '■) + ^« + - + "I - '°g(i/.)) 

logf , "" + '■)(' + '') \ (6,20) 



^{t + a){t + b) \2^{t + a){t + b) + 2t + a + b^ 

vr— , t^oo. (6.21) 



Note that the correct large t behavior is reproduced in (6.21). Some trivial algebra yields (6.2). 

Now consider (6.6)-(6.9). We will explicitly derive (6.9); the integral (6.8) is then obtained 

by the analytical continuation of (6.9) to t = —1, whereas the integrals (6.6) and (6.7) are 

obtained by taking t — )■ oo and t — )■ in (6.9) respectively. With the Schwinger parametrization 

70 



(6.15) and the partial fraction decomposition 

X 1/A 



+ 



(x-l/A)(x + t) (x-A)(t + l/A) {t + l/\){x + ty 
the integral becomes 

log(l -x) J _ [^ dX /-^ / 1/A t \ dx 



a {x + t)^{b-x){x-a) Jo l + ^tj^ \x-l/\ x + tj ^(^b-x)ix-a) 



"V 1 t \ dX 

vr / = + 



V y'il-XaJil-Xb) ^{t + a){t + h)) 1 + At' 

The last equation was obtained by invoking (6.16) and taking the analytic continuation of t to 
— 1/A, together with the implicit assumption that 3f? t > in (6.16). From a further change of 
variable 1 + A t = x, we have 

log(l - x) , /•^+* ( I 1 \ dx 

dx = TT 



a {x + t)y^{b-x){x-a) Ji \^y{t + a){t + b) y'{t + a- ax){t + b - bx) ) x 

log(l + t) 



vr 



^{t + a){t + b) 



2^J{t + a){t + b)^{t + a-ax){t + b-bx)- [t{a + b)+ 2ab]x + 2(t + a){t + b) 



x=l+t-[ 



x=l 



vr , 2J(t + a)(t + b)J{l - a){l - b) + \2 - a - b]t + a + b -2ab 

log- 



v/(t + a)(i + 6) 2^J{t + a){t + b) + 2t + a + b 

vr, 2J{l-a){l-b) + 2-a-b 
~-log ^ , t-^oo. 

Simple algebra yields (6.9). 

Finally, consider the principal value integrals (6.10) and (6.11). These results are obtained 
by taking the square root to the denominator of the integrand, followed by performing a partial 
fraction decomposition and invoking (6.18) and (6.19). 
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